#include <stdio.h>
#include <time.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define MAX_POINTS      1000000
#define MAX_DISTANCE    1000000
//#define K1              0               //T1 parameter
//#define K2              0               //T2 parameter
//#define RAND_MAX        32767
//#define RURAL           10              //Parameter used to control the randomrate of rural areas.

/*Global Variables*/
float K1,K2;
float A1,A2;
float B1,B2;
int RURAL;


unsigned short consumerNumber;
float branchRate;
int ViSize;
float tmpX;
float tmpY;
int fracalCount;
float NetworkLength;

/*Claim functions*/
int checkInput(int cousumerNumber, float branchRate);
void getVerticalPoint(float x1,float y1,float x2,float y2,float x3,float y3,float *positionX,float *positionY,int *sgn,float *distance);
float getNetworkLength(float* A,int N,int consumerNumber);
void rearrangePoints(float V[][2]);
int generateConnection(float A[][2*consumerNumber],float V[][2],float Vi[][2]);
void wirte2File(float A[][2*consumerNumber],float V[][2],float Vi[][2]);
void generateConsumerPosition(float A[][2*consumerNumber],float V[][2],float Vi[][2]);
int getClosestPoint(float V[][2],int num);
int getLinkedPoint(float A[][2*consumerNumber],int num,int selfPoints);

/*Main Entrance*/
int main(int argc, char* argv[])
{
    //Read input numbers
    consumerNumber = atof(argv[1]);
    branchRate = atof(argv[2]);
    K1 = atof(argv[3]);
    K2 = atof(argv[4]);
    A1 = atof(argv[5]);
    A2 = atof(argv[6]);
    RURAL = atoi(argv[7]);
    B1 = atof(argv[8]);
    B2 = atof(argv[9]);
    if  ( checkInput(consumerNumber,branchRate)==1  ) return 0;

    //Init...
    float V[consumerNumber][2];                     //Matrix used to store consumer positions
    float A[2*consumerNumber][2*consumerNumber];	//Linking Matrix
    float Vi[2*consumerNumber][2];                  //Positions including T-joints

    //Generate Consumer positions
    generateConsumerPosition(A,V,Vi);


    //Rearrange points for different Branchrate
    rearrangePoints(V);

    //Generate a simple connection
    generateConnection(A,V,Vi);



    //print the Networklength
    NetworkLength = getNetworkLength(A,ViSize,consumerNumber);  //Change the definition of it.
    printf("\nNetworkLength:%f\n",NetworkLength);

    //export the generated Linking Matrix to A.txt , V.txt, Vi.txt
    wirte2File(A,V,Vi);
    return 0;
}

int checkInput(int consumerNumber, float branchRate)
{
    if (consumerNumber>10 && consumerNumber<MAX_POINTS && branchRate>0 && branchRate<1)
    {
        return 0;
    }
    else
    {
        printf("%s\n", strerror(5));	//Input Error
        return 1;
    }
}


void getVerticalPoint(float x1,float y1,float x2,float y2,float x3,float y3,float *positionX,float *positionY,int *sgn,float *distance)
{

    if ((x3-x1)*(x2-x1)+(y3-y1)*(y2-y1)>0 && (x3-x2)*(x1-x2)+(y3-y2)*(y1-y2)>0)
    {
        *sgn=1;
    }
    else
    {
        *sgn=0;
    }
    //Ax+By+C=0
    float A =  (y2-y1)/(x2-x1);
    float B = -1;
    float C =  y1 - A*x1;
    float tmp = A*A+B*B;
    *positionX =	(B*B*x3-A*B*y3-A*C)/tmp;
    *positionY = 	(A*A*y3-A*B*x3-B*C)/tmp;
    *distance = (*positionX - x3)*(*positionX - x3) + (*positionY - y3)*(*positionY - y3);


}


float getNetworkLength(float* A,int N,int consumerNumber)
{

    float D[N][N];
    int i,j,k;
    double Asum=0;
    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            if ((*(A+i*2*consumerNumber+j))==0)
                D[i][j]=99999;
            else
                D[i][j]=*(A+i*2*consumerNumber+j);

        }
    }


     for(i=0;i<N;i++)
     {
         D[i][i]=0;
     }

     for (k=0;k<N;k++)            //Floyd
     {

         for (i=0;i<N;i++)
         {
             for (j=0;j<N;j++)
             {
                 if (D[i][j]> D[i][k]+D[k][j])
                    D[i][j]=D[i][k]+D[k][j];
             }
         }
     }


     for(i=0;i<N;i++)
     {
        for(j=0;j<N;j++)
        {
           Asum = Asum+D[i][j];
        }
     }

     float aver_D = (float) (Asum/(N*(N-1)));
     if (aver_D>999)
     {
         printf("Not a linked gprah");
         return -1;
     }
     else
         return aver_D;

}


//Rearrange points for different Branchrate
void rearrangePoints(float V[][2])
{
    float length;
    int num;
    int i,j;
    float tmp[2];
    for (i=2;i<(int)(branchRate*consumerNumber)+1;i++)
    {
    length = MAX_POINTS;
    num = 0;
    for (j=i;j<consumerNumber+1;j++)
    {
        if (	((V[j-1][0]-V[i-2][0])*(V[j-1][0]-V[i-2][0])+(V[j-1][1]-V[i-2][1])*(V[j-1][1]-V[i-2][1]) ) <length	)
        {
            length = ((V[j-1][0]-V[i-2][0])*(V[j-1][0]-V[i-2][0])+(V[j-1][1]-V[i-2][1])*(V[j-1][1]-V[i-2][1]) ) ;
            num = j;
        }
    }
        tmp[0] = V[i-1][0];
        tmp[1] = V[i-1][1];
        V[i-1][0]=V[num-1][0];
        V[i-1][1]=V[num-1][1];
        V[num-1][0]=tmp[0];
        V[num-1][1]=tmp[1];
    }
    printf("%d consumer points generated successfully, with a branchRate of%0.2f.\n",consumerNumber,branchRate);
}

//Generate a SimpleConnection
int generateConnection(float A[][2*consumerNumber],float V[][2],float Vi[][2])
{
    int i,j;
    int inputSize = consumerNumber;
    for(i=0;i<2*consumerNumber;i++)					//Init A- The linking Matrix
    {
        for(j=0;j<2*consumerNumber;j++)
        {
            A[i][j]=0;
        }
    }
    //Init Vi- The points Matrix containing T-points.
    for(i=0;i<consumerNumber;i++)
    {
        Vi[i][0]=V[i][0];
        Vi[i][1]=V[i][1];
    }
    int a1,a2,a1Ref,a2Ref,sgn;
    float positionX,positionY,positionXRef,positionYRef,distance,distanceRef;
    ViSize = consumerNumber;

    //Generating
    for(i=2;i<inputSize+1;i++)
    {
        if(i==2)
        {
            A[0][1]=1;
            A[1][0]=1;
        }
        else
        {
            distanceRef=MAX_DISTANCE;
            for(a1=1;a1<i;a1++)
            {
                for(a2=a1+1;a2<i;a2++)
                {
                    //Get VerticalPoint
                    getVerticalPoint(V[a1-1][0],V[a1-1][1],V[a2-1][0],V[a2-1][1],V[i-1][0],V[i-1][1],&positionX,&positionY,&sgn,&distance);
                    if (sgn==1 && distance<distanceRef && A[a1-1][a2-1]==1)
                    {
                        positionXRef=positionX;
                        positionYRef=positionY;
                        distanceRef=distance;
                        a1Ref=a1;
                        a2Ref=a2;
                    }
                }
            }
                for(a1=1;a1<i;a1++)
                {
                    if (	(Vi[i-1][0]-Vi[a1-1][0])*(Vi[i-1][0]-Vi[a1-1][0])+(Vi[i-1][1]-Vi[a1-1][1])*(Vi[i-1][1]-Vi[a1-1][1]) < distanceRef )
                    {
                        positionXRef=Vi[a1-1][0];
                        positionYRef=Vi[a1-1][1];
                        distanceRef=(Vi[i-1][0]-Vi[a1-1][0])*(Vi[i-1][0]-Vi[a1-1][0])+(Vi[i-1][1]-Vi[a1-1][1])*(Vi[i-1][1]-Vi[a1-1][1]);
                        a1Ref=a1;
                        a2Ref=a1;
                    }
                }
                for(a1=1;a1<i;a1++)
                {
                    for (a2=inputSize+1;a2<ViSize+1;a2++)
                    {
                        getVerticalPoint(Vi[a1-1][0],Vi[a1-1][1],Vi[a2-1][0],Vi[a2-1][1],Vi[i-1][0],Vi[i-1][1],&positionX,&positionY,&sgn,&distance);
                        if (sgn==1 && distance<distanceRef && A[a1-1][a2-1]==1)
                        {
                            positionXRef=positionX;
                            positionYRef=positionY;
                            distanceRef=distance;
                            a1Ref=a1;
                            a2Ref=a2;
                        }
                    }
                }

                if (a1Ref==a2Ref)
                {
                    A[i-1][a1Ref-1]=1;
                    A[a1Ref-1][i-1]=1;
                }
                else if ((V[i-1][0]-positionXRef)*(V[i-1][0]-positionXRef)+(V[i-1][1]-positionYRef)*(V[i-1][1]-positionYRef)<0.00001)
                {

                    A[a1Ref-1][i-1]=1;
                    A[i-1][a1Ref-1]=1;
                    A[a2Ref-1][i-1]=1;
                    A[i-1][a2Ref-1]=1;
                    A[a2Ref-1][a1Ref-1]=0;
                    A[a1Ref-1][a2Ref-1]=0;
                }
                else
                {
                    ViSize++;
                    Vi[ViSize-1][0] = positionXRef;
                    Vi[ViSize-1][1] = positionYRef;
                    A[i-1][ViSize-1]=1;
                    A[ViSize-1][i-1]=1;
                    A[a1Ref-1][ViSize-1]=1;
                    A[ViSize-1][a1Ref-1]=1;
                    A[a2Ref-1][ViSize-1]=1;
                    A[ViSize-1][a2Ref-1]=1;
                    //Disconnect The link
                }
        }
    }

    return ViSize;
}


void wirte2File(float A[][2*consumerNumber],float V[][2],float Vi[][2])
{
    int i,j;
    FILE *fp;
    fp=fopen("A.txt","w");
    for(i=0;i<ViSize;i++)
    {
        for(j=0;j<ViSize;j++)
        {
            fprintf(fp,"%f,",A[i][j]);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);

    FILE *fp1;
    fp1=fopen("Vi.txt","w");
    for(i=0;i<ViSize;i++)
    {
        fprintf(fp1,"%f,%f\n",Vi[i][0],Vi[i][1]);

    }
    fclose(fp1);

    FILE *fp2;
    fp2=fopen("V.txt","w");
    for(i=0;i<consumerNumber;i++)
    {
        fprintf(fp2,"%f,%f\n",V[i][0],V[i][1]);

    }
    fclose(fp2);

    FILE *fp3;
    fp3=fopen("F.txt","w");
    fprintf(fp3,"%f", NetworkLength);
    fclose(fp3);
}


//Generate consumer position. based on the given t1,t2,branchrate,and number.
void generateConsumerPosition(float A[][2*consumerNumber],float V[][2],float Vi[][2])
{
    int i;
    double t1,t2;;
    int closestRow,linkedRow;
    V[0][0]=0.5;
    V[0][1]=0.5;
    srand((unsigned)time(0));
    fracalCount=0;
    for (i=1;i<4;i++)
    {
        V[i][0]=(float)( rand())/RAND_MAX;
        V[i][1]=(float)( rand())/RAND_MAX;
    }
    for (i=4;i<consumerNumber;i++)
    {
        t1=K1*pow((double)(i+1),A1)-B1*(i+1)/20;
        t2=K2*pow((double)(i+1),A2)+B2*(i+1)/20;
        V[i][0]=(float)( rand())/RAND_MAX;
        V[i][1]=(float)( rand())/RAND_MAX;
        ViSize = generateConnection(A,V,Vi);
        float len = getNetworkLength(A,ViSize,consumerNumber);
        if(len>t2 || len<t1)
        {
            //Find a repalcement point
            closestRow = getClosestPoint(V,i);
            linkedRow = getLinkedPoint(A,closestRow,i);
            //if (closestRow ==0)
            //    closestRow =1;

            float positionX,positionY,distance;
            int sign;
            //getVerticalPoint(V[closestRow][0],V[closestRow][1],V[closestRow-1][0],V[closestRow-1][1],V[i][0],V[i][1],&positionX,&positionY,&sign,&distance);
            getVerticalPoint(V[closestRow][0],V[closestRow][1],V[linkedRow][0],V[linkedRow][1],V[i][0],V[i][1],&positionX,&positionY,&sign,&distance);

            //Template
            //printf("%d,%d, %f, %f, %f, %f, %f, %f\n",i,closestRow,positionX,positionY,V[i][0],V[i][1],V[closestRow][0],V[closestRow][1]);
            fracalCount++;
            //Make sure that the point is in side the border
            if(positionX<0.01 || positionY<0.01 ||positionX>1 ||positionY>1 )
            {
                V[i][0]=(float)( rand())/RAND_MAX;
                V[i][1]=(float)( rand())/RAND_MAX;

            }
            //Make sure that there are no same points(which can lead to a big problem)
            else if(positionX==V[closestRow][0] && positionY==V[closestRow][1])
            {
                V[i][0]=(float)( rand())/RAND_MAX;
                V[i][1]=(float)( rand())/RAND_MAX;
            }
            else
            {
                V[i][0]=positionX + (rand())/RAND_MAX/RURAL-1/RURAL/2;
                V[i][1]=positionY + (rand())/RAND_MAX/RURAL-1/RURAL/2;
            }
        }
    }
    printf("FracalCount:%d\n",fracalCount);

}

int getClosestPoint(float V[][2],int num)
{
    int i;
    float length = MAX_DISTANCE;
    int point = 0;
    float x = V[num][0];
    float y = V[num][1];
    for (i=0;i<num-1;i++)
    {
        if ((x-V[i][0])*(x-V[i][0]) + (y-V[i][1])*(y-V[i][1]) < length)
        {
            length = (x-V[i][0])*(x-V[i][0]) + (y-V[i][1])*(y-V[i][1]) ;
            point = i;
        }
    }
    return point;
}

int getLinkedPoint(float A[][2*consumerNumber],int num, int selfPoints)
{
    int i;
    for (i=num-1;i>0;i=i-1)
    {
        if (A[num][i]==1)
        {
            return i;
        }
    }
    return num-1;
}
